﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using MathNet.Numerics.LinearAlgebra;
using MentalAlchemy.Atomics;

namespace MentalAlchemy.Molecules.MachineLearning
{
	/// <summary>
	/// [molecule]
	/// 
	/// AMaLGaM IDEA implementation with anticipated mean shift (AMS) using description from:
	///		Bosman, P.; Grahl, J. & Thierens, D. (2007), 
	///		'Adapted Maximum-Likelihood Gaussian Models for Numerical Optimization with Continuous EDAs'
	///		Technical report, CWI
	/// see p. 14. The lines numbers are set according to that algorithm's pseudocode.
	/// </summary>
	public class AmalgamIDEA : EvolutionaryAlgorithm
	{
		protected const float tau = 0.3f;	// truncation selection threshold.
		protected float thetaSdr = 1;	// threshold for normalized improvements.
		protected const float nDecrease = 0.9f;	// mean shift multiplier
		protected const float nIncrease = 1 / nDecrease;	// mean shift multiplier
		protected float alpha = tau / (2 - 2*tau);
		protected float delta = 2;
		protected float cMult;
		protected int t;
		float[] mean, muPrev;
		float[,] cov;

		protected List<Individual> offsp;

		public override void Init(EAParameters parameters)
		{
			parameters.MinGeneValue = -10;
			parameters.GeneValueRange = 15;
			base.Init(parameters);
			cMult = 1;
		}

		public override void Continue(EAParameters parameters)
		{
			//
			// validate [parameters].
			if (FitnessFunction != null && !EAElements.ValidateParameters(parameters)) throw new Exception("[AmalgamIDEA.Continue]: Invalid parameters setting or fitness function is undefined.");

			Evaluate();

			//
			// Main cycle.
			//var muPrev = new float[1];
			int selSize = (int)(tau*popul.Count);
				//,
				//offspSize = popul.Count - selSize,
				//shiftSize = (int)(alpha * selSize);
			for (CurrentGeneration = 1; CurrentGeneration <= parameters.GenerationsNumber; ++CurrentGeneration)
			{
				//selPopul = EAElements.TruncationSelection(popul, selSize);	// line 11.01
				Select(parameters);

				// compute ML estimates.
				//float[] mean;
				//float[,] cov;
				//EAElements.ComputeMLEstimates(selPopul, out mean, out cov);	// line 11.02
				//cov = MatrixMath.Mul(cov, cMult);	// line 11.03
				//var gaussDistr = new GaussianDistribution(mean, cov);

				////
				//// create offsprings.
				//var offsp = new List<Individual>();
				//for (int i = 0; i < offspSize; i++)	// 11.04
				//{
				//    var tempOffsp = new Individual();
				//    tempOffsp.Genes = new List<float>(gaussDistr.Next());
				//    offsp.Add(tempOffsp);	// 11.04.1
				//}
				Cross(parameters);

				//
				// apply anticipated mean shift.
				//if (CurrentGeneration > 1)
				//{
				//    var muShift = VectorMath.Sub(mean, muPrev);	//	11.05.1
				//    var shift = VectorMath.Mul(muShift, cMult * delta);
				//    for (int i = 0; i < shiftSize; i++)
				//    {
				//        offsp[i].Genes = new List<float>(VectorMath.Add(offsp[i].Genes.ToArray(), shift));	// 11.05.2.1
				//    }
				//}

				////
				//// find current best.
				//var bestFit = EAElements.GetBestIndividual(selPopul).Fitness;	// 11.06

				////
				//// Find out about improvements by offspring.
				//var nImprovement = 0;	// 11.07. Number of cases when offspring is better than parents.
				//var xImprovement = VectorMath.Zeros(offsp[0].Genes.Count);	// 11.08	// average direction of improvement by offsprings.
				//for (int i = 0; i < offspSize; i++)
				//{
				//    offsp[i].Fitness = FitnessFunction.Compute(offsp[i].Genes.ToArray());	// 11.09.1
				//    if (FitnessComparator.IsBetter(offsp[i].Fitness, bestFit))	// 11.09.2
				//    {
				//        nImprovement++;
				//        VectorMath.Accumulate(ref xImprovement, offsp[i].Genes.ToArray());
				//    }
				//}

				////
				//// apply offpsring improvement.
				//if (nImprovement > 0)	// 11.10
				//{
				//    VectorMath.Mul(ref xImprovement, 1f/nImprovement);	// 11.10.1

				//    //// get Cholesky decomposition matrix.
				//    //var m = Matrix.Create(MatrixMath.ConvertToDoubles(cov));
				//    //var cholTri = m.CholeskyDecomposition.TriangularFactor;	// todo: check this!
				//    //var inv = cholTri.Inverse();	// invert Cholesky triangular factor.

				//    // compute rho
				//    var rho = 0f;
				//    for (int i = 0; i < parameters.IndividualSize; i++)	// 11.10.2
				//    {
				//        var temp = Math.Abs(xImprovement[i] - mean[i])/cov[i, i];
				//        if (temp > rho) { rho = temp; }
				//    }

				//    if (rho > thetaSdr)
				//    {
				//        cMult *= nIncrease;	// 11.10.3.1
				//    }
				//}
				//else
				//{
				//    cMult *= nDecrease;	// 11.10.1
				//}

				//if (cMult < 1){cMult = 1;} // 11.11 & 11.11.2
				Mutate(parameters);

				//popul.Clear();
				//popul.AddRange(selPopul);
				//for (int i = 0; i < offsp.Count; i++) {popul.Add(offsp[i]);}	// 11.12 & 11.13 because fitness is copied with genes.

				//muPrev = mean;	// 11.14

				//CollectStats();
				NextGeneration(parameters);
			}
		}

		/// <summary>
		/// Selection.
		/// </summary>
		/// <param name="parameters"></param>
		public override void Select(EAParameters parameters)
		{
			int selSize = (int)(tau * popul.Count);
			selPopul = EAElements.TruncationSelection(popul, selSize);	// line 11.01
		}

		/// <summary>
		/// Crossing
		/// </summary>
		/// <param name="parameters"></param>
		public override void Cross(EAParameters parameters)
		{
			var offspSize = popul.Count - selPopul.Count;

			// compute ML estimates.
			EAElements.ComputeMLEstimates(selPopul, out mean, out cov);	// line 11.02
			cov = MatrixMath.Mul(cov, cMult);	// line 11.03
			var gaussDistr = new GaussianDistribution(mean, cov);

			//
			// create offsprings.
			offsp = new List<Individual>();
			for (int i = 0; i < offspSize; i++)	// 11.04
			{
				var tempOffsp = new Individual();
				tempOffsp.Genes = new List<float>(gaussDistr.Next());
				offsp.Add(tempOffsp);	// 11.04.1
			}
		}

		/// <summary>
		/// Mutation and update of multipliers.
		/// </summary>
		/// <param name="parameters"></param>
		public override void Mutate(EAParameters parameters)
		{
			var shiftSize = (int)(alpha * selPopul.Count);
			var offspSize = popul.Count - selPopul.Count;

			//
			// apply anticipated mean shift.
			if (CurrentGeneration > 1)
			{
				var muShift = VectorMath.Sub(mean, muPrev);	//	11.05.1
				var shift = VectorMath.Mul(muShift, cMult * delta);
				for (int i = 0; i < shiftSize; i++)
				{
					offsp[i].Genes = new List<float>(VectorMath.Add(offsp[i].Genes.ToArray(), shift));	// 11.05.2.1
				}
			}

			//
			// find current best.
			var bestFit = EAElements.GetBestIndividual(selPopul).Fitness;	// 11.06

			//
			// Find out about improvements by offspring.
			var nImprovement = 0;	// 11.07. Number of cases when offspring is better than parents.
			var xImprovement = VectorMath.Zeros(offsp[0].Genes.Count);	// 11.08	// average direction of improvement by offsprings.
            for (int i = 0; i < offspSize; i++)
			{
				offsp[i].Fitness = FitnessFunction.Compute(offsp[i].Genes.ToArray());	// 11.09.1
				if (FitnessComparator.IsBetter(offsp[i].Fitness, bestFit))	// 11.09.2
				{
					nImprovement++;
					VectorMath.Accumulate(ref xImprovement, offsp[i].Genes.ToArray());
				}
			}

			//
			// apply offpsring improvement.
			if (nImprovement > 0)	// 11.10
			{
				VectorMath.Mul(ref xImprovement, 1f/nImprovement);	// 11.10.1

				//// get Cholesky decomposition matrix.
				//var m = Matrix.Create(MatrixMath.ConvertToDoubles(cov));
				//var cholTri = m.CholeskyDecomposition.TriangularFactor;	// todo: check this!
				//var inv = cholTri.Inverse();	// invert Cholesky triangular factor.

				// compute rho
				var rho = 0f;
				for (int i = 0; i < parameters.IndividualSize; i++)	// 11.10.2
				{
					var temp = Math.Abs(xImprovement[i] - mean[i])/cov[i, i];
					if (temp > rho) { rho = temp; }
				}

				if (rho > thetaSdr)
				{
					cMult *= nIncrease;	// 11.10.3.1
				}
			}
			else
			{
				cMult *= nDecrease;	// 11.10.1
			}

			if (cMult < 1){cMult = 1;} // 11.11 & 11.11.2
		}

		/// <summary>
		/// Advance to the next generation.
		/// </summary>
		/// <param name="parameters"></param>
		public override void NextGeneration(EAParameters parameters)
		{
			popul.Clear();
			popul.AddRange(selPopul);
			for (int i = 0; i < offsp.Count; i++) { popul.Add(offsp[i]); }	// 11.12 & 11.13 because fitness is copied with genes.

			muPrev = mean;	// 11.14

			CollectStats();
		}

		///// <summary>
		///// 
		///// This is a C# version of Peter Bosman's funciton:
		/////		double *generateNewSolution( int population_index )
		///// from the file "AMaLGaM-Bayesian.c"
		///// </summary>
		///// <returns></returns>
		//public static float[] GenerateNewSolution (int paramsCount, int popIndex)
		//{
		//    int i, j, var_index, times_not_in_bounds;
		//    float[] result = new float[paramsCount];
		//    float mean;

		//    //(double *) Malloc( number_of_parameters*sizeof( double ) );

		//    for( i = 0; i < paramsCount; i++ )
		//    {
		//        var_index = sampling_ordering[popIndex][i];

		//        mean = mean_vectors[population_index][var_index]*precision_matrices_first_col[population_index][var_index][0];
		//    for( j = 0; j < parent_vectors_lengths[population_index][var_index]; j++ )
		//      mean -= (result[parent_vectors[population_index][var_index][j]] - mean_vectors[population_index][parent_vectors[population_index][var_index][j]])*precision_matrices_first_col[population_index][var_index][j+1];

		//    if( precision_matrices_first_col[population_index][var_index][0] == 0 )
		//      mean = mean_vectors[population_index][var_index];
		//    else
		//      mean /= precision_matrices_first_col[population_index][var_index][0];

		//    if( variance_vectors_sampling[population_index][var_index] <= 0.0 )
		//    {
		//      result[var_index] = mean;
		//      if( !isParameterInRangeBounds( result[var_index], var_index ) )
		//        result[var_index] = lower_init_ranges[var_index] + (upper_init_ranges[var_index] - lower_init_ranges[var_index])*randomRealUniform01();
		//    }
		//    else
		//    {
		//      times_not_in_bounds = -1;
		//      out_of_bounds_draws[population_index]--;
			  
		//      do
		//      {
		//        times_not_in_bounds++;
		//        samples_drawn_from_normal[population_index]++;
		//        out_of_bounds_draws[population_index]++;
			    
		//        if( times_not_in_bounds >= 100 )
		//          result[var_index] = lower_init_ranges[var_index] + (upper_init_ranges[var_index] - lower_init_ranges[var_index])*randomRealUniform01();
		//        else
		//          result[var_index] = random1DNormalParameterized( mean, variance_vectors_sampling[population_index][var_index] );
		//      }
		//      while( !isParameterInRangeBounds( result[var_index], var_index ) );
		//    }
		//    }

		//    return( result );
		//}
	}
}
